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Abstract 



This paper focuses on large neural networks whose synaptic connectivity matrices are 
randomly chosen from certain random matrix ensembles. The dynamics of these networks 
can be characterized by the eigenvalue spectra of their connectivity matrices. In reality, 
neurons in a network do not necessarily behave in a similar way, but may belong to sev- 
eral different categories. The first study of the spectra of two-component neural networks 
was carried out by Rajan and Abbott. In their model, neurons are either 'excitatory' or 
'inhibitory', and strengths of synapses from different types of neurons have Gaussian dis- 
tributions with different means and variances. A surprising finding by Rajan and Abbott 
is that the eigenvalue spectra of these types of random synaptic matrices do not depend 
on the mean values of their elements. In this paper we prove that this is true even for a 
much more general type of random neural network, where there is a finite number of types 
of neurons, and their synaptic strengths have correlated distributions. Furthermore, using 
the diagrammatic techniques, we calculate the explicit formula for the spectra of synaptic 
matrices of multi-component neural networks. 

1 Introduction 

In neuroscience, interconnections of neurons are often represented by synaptic matrices whose 
elements are drawn from a certain random matrix ensemble [TJ [5] . Knowing the distribution 
of eigenvalues of these random matrices is very important in studying spontaneous activities 
and evoked responses of the network. To calculate the eigenvalue distribution of these ma- 
trices, it is often necessary to work with asymmetric (non-hermitean) random matrix theory, 
which has been successfully applied to many fields of physics and interdisciplinary sciences, 
e.g. the phase diagram of QCD [3j 2] , nuclear decay and resonances in multichannel chaotic 
scattering [5] and neural networks [TJ EJ |6j. 

A prominent result of asymmetric random matrix theory is Girko's circle law [7]. In 
its variation with partial symmetry, the circle becomes an ellipse [5J. These classic results, 
however, can not be directly applied to realistic neural network models where neurons do 
not behave in the same way [TJ [8J [9] . Assume there are N number of neurons and let W 
be the synaptic matrix. In the model of Rajan and Abbott [JJ, there are fN number of 
neurons which are 'excitatory', and all others are 'inhibitory'. To model this neural network, 



elements in fN columns of W are sampled from a Gaussian distribution with mean /j,e/\/N 
and variance <J E /N and elements in the remaining (1 — f)N columns of W are Gaussian 
variables with mean /j,i/y/N and variance crj/N. Therefore, the synaptic matrix has the 
structure W = J A + M, where J is drawn from the real Ginibre ensemble 12 such that 
(Jij)=0, (Jfj) — l/N, and A = di&g(<jElfN, cr//(i_/)jv)i where IfN and I^_f^ N are identity 
matrices of dimension fN and (l — f)N, respectively. M is a constant matrix whose elements 
are the mean strength of the synapses. Since there are two types of synapses, every row of M 
is identical and in each row, the first fN elements are equal to (J,e/V~N~ and the remaining 
(1 — f)N elements equal to /i/ / y/N. In particular, M is chosen to be in a 'balanced' situation 
such that fuE + (1 — /)/«/ = [131 E]. To confine eigenvalues inside a unit circle, a second 
constraint [T] is introduced which requires that the strengths of the synapses attached to each 
neuron independently sum to zero. It is found in [1] that, in the limit N — > oo, modifying the 
mean strengths of excitatory and inhibitory synapses has no effect on the eigenvalue spectra 
of the synaptic matrices. Therefore, the spectrum of W is identical to that of J A. 

It is natural to wonder why the 'mean' strength matrix M has no effect on the spectra. 
Moreover, in real biological neural systems, several different types of neurons may connect 
each other to form a multi-component network [8l [9] . Distributions of synaptic strengths of 
different types of neurons are distinct [10] and could be non-Gaussian [TIB [11] . Dynamics of 
this network therefore depend on properties of each type of neuron. It is interesting to find 
whether or not the eigenvalue density of this type of networks depends on the mean value 
of each individual type of synapse. These questions are addressed in section 2. One of the 
main results of this paper shows that even without the second constraint in |T| and when 
synaptic strengths have certain non-Gaussian distributions, the spectrum of the network still 
does not depend on the mean synaptic strengths. 

Aside from biological motivations, the eigenvalue problem of random plus fixed matrices 
has been a research topic in both random matrix theory and condensed matter physics for a 
long time [TSJ [TSJ [T71 [T5]. A different point of view of the problem in this paper is: how is 
the density function of large random matrix J A perturbed by the rank-1 constant matrix M. 
Note that random matrix J A is not of Wigner type, nor are its elements independently and 
identically distributed (iid). Therefore this paper provides new results to similar problems 
studied in [I7| US]. 

Furthermore, finding the eigenvalue density of random matrices of the form J A, where 
J is drawn from a random matrix ensemble and A is a fixed matrix, has been an interesting 
topic in random matrix theory and mesoscopic physics |19j . When J is drawn from the 
circular unitary ensemble (CUE), an exact result is given in [20], and the large-N limit is 
calculated in [21]. In section 3, we calculate the density function of J A where J belongs to 
the real Ginibre ensemble using the method introduced in [221 123) . Discussion and remarks 
are made in the last section. 



2 Synaptic strength of non-Gaussian distributions 

Let the N x N dimensional real matrix W be the synaptic matrix of an A^-neuron network. 
Assume there are m types of neurons and the i-th type of neuron has a population of fiN, 
Y^iLi fi = 1- Define a constant diagonal matrix 

A = diag(cri// 1 jv, • ■ • , (T m I fmN ) > 0, (2.1) 

where Ifcjsr is the /^Af-dimensional identity matrix. Let v be an A^-dimensional row vector 
with the following form, 

V = ( ^1, ■ ■ ■ , Vl , V2, ■ ■ ■ , ■ ■ ■ , Vm, ■ ■ ■ , Vm) , (2-2) 
fiN f 3 N f m N 
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where [ii is the mean strength of the synapses from neurons of the i-th type. Define the 
N x N dimensional 'mean' matrix M, whose rows are all equal to v. The synaptic matrix 
W in our model takes the form 

W = J A + M, (2.3) 
where J is an N x N dimensional real random matrix drawn from the ensemble 

P(J) = | exp(-iVtrT/(JJ T )), (2.4) 

where V is an arbitrary function and Z is the normalization constant. The case V(x) = x/2 
corresponds to Ginibre ensemble where elements of J are statistically independent Gaus- 
sian variables. By symmetry, the mean of Jij vanishes, (Jij) = 0. The variance of Jy is 
determined by V, which is normalized to be 1/N, i.e. (Jfj) = 1/N. 

Synaptic matrix W defined in Eq. (|2.3|> has m column blocks. Each block corresponds to 
one type of neuron. By construction, elements of W have the following statistical properties, 

Var(W^) = Ajj/N and (W zj ) = v v (2.5) 

i.e. the i-th type of synaptic strength has variance crf/N and mean fa. The matrix M has 
a similar column block structure as the one defined in [T]. Following p], we also choose to 
put the synapses at the 'balanced' situation, i.e. 

m 

£/iA*i = 0. (2.6) 

We want to know the eigenvalue density p(x,y) of W in the limit N — > oo, with fixed 
/'s. Let z = x + iy, the density p{x, y) is related to the Green's function Gw{z, z) as 

P{x,y) = --^G w (z,z), where G w (z,z) = i (tr N ) . (2.7) 

7T oz N \ z — W / j 

In the above formula, (• • • ) j means averaging over the ensemble Eq. ([2.4[) of matrix J. We 
write the Green's function as Gw(z, z) to emphasize it is not analytic on a 2-dimensional 
region of the (x, y)-plane, more details can be found in [23 . This region is called the support 
of the density function since on which we have J^Gw{z, z) ^= 0. Since we will be dealing 
with both N x N and 2iV x 2N dimensional matrices, to remove ambiguity, we use trjy as 
the trace operator for N x TV matrices. We will work on asymmetric random matrices with 
the methods introduced in |22[ [23] . For consistency, we adopt the notation convention of 
[22] in the remaining of this paper. Define a 2N x 2N dimensional matrix 

Z=(l iV (2-8) 



z 



For asymmetric matrix W, define the resolvent (matrix valued Green's function) Qw 

as 



Qw{Z) 



Gi Qi 



W 



W 



T 



(2.9) 



Introducing the self-energy Ey^, we have 



Qw{Z) = - l - ■ (2.10) 



Z-T, 
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The Green's function in Eq. (|2.7|) can be found from Gw [22l I23j . 



G w (z,z)= lim — trjv£?i, (2-11) 

where the limit N — > oo is taken before A — > 0. Similarly, as in Eqs. (|2.7|) - (12.1ip , we define 
Gja, Sja and Gja- Introducing a constant matrix 

M=( M MT ), (2.12) 



we have the relation 



Qw{Z) = Gja(Z -M)= ! . (2.13) 

j6 — A4 — 2jja 

It is impossible to calculate Eja explicitly for arbitrary V. But for our purpose it is sufficient 
to know its basic structure. Without loss of generality, assume V{x) = x/2 + . . . , so that we 
can expand tiV( JJ T ) as 

tr N V(JJ T ) = ^ti N JJ T + g 2 tT N (JJ T ) 2 + g 3 tT N {JJ T ) 3 + ■■■ . (2.14) 

We expand the higher order terms in Eq. (|2.14j ) and use the quadratic term to calculate the 
ensemble averages, denoted by ( )o- Let W = JA. Then because of the presence of matrix 
A, we have 

(W ab W? d ) = ^S ad 6 bc , (W ab ) =0. (2.15) 




Figure 1: Contributions of the quadratic cumulant T2 and quartic cumulant T4 to the self- 
engergy S JA . 

The self-energy Sja can be written in terms of the cumulants of P(J), i.e. I^/c, k — 
1,2,.... And it is well known that in the limit N — > 00, to leading order in i, each of these 
cumulants is the sum of all connected plannar diagrams with k external J's and J T 's. All 
diagrams which contribute to Eja are also planar diagrams, as shown in Figll] By Eq. (|2.15|) . 
the self-energy Eja has the following structure 

where scalars a and b are functions of z and z and are determined by V and A is defined 
in Eq. (|2.ip . In appendix A, we prove that Gw = Gja- This fact, together with Eq. (|2.7l) . 
completes the proof that the eigenvalue spectrum of W is identical to that of the random 
matrix JA, as discovered in IT] when J belongs to Ginibre ensemble. 



4 



300 



^ 200 
CO 

c 

O 100 



300 



300 



£,200 
CO 

c 

^ 100F 



0.2 
|z| 



0.2 
|Z| 




Figure 2: Density p of eigenvalues as a function of radius \z\ in the complex plane. Sim- 
ulations are run for random matrices of dimension N = 200, drawn from the ensemble de- 
fined in Eq. (|2.4p with V{x) = x + x 2 . There are four types of neurons in the network, with 
/ = (0.1,0.2,0.3,0.4) and A = diag(0.5J 2 o, 1-0/40, 1.5/bo, 2.0J 80 ). Panel 1, p = (10,30,30,-40). 
Panel 2, fj, = (1,3,3,-4). Panel 3, p = (0,0,0,0), i.e. M = 0. Panel 4, density functions in 
panel 1-3 are drawn in the same plane. We find the bulk of the three functions are very similar. 
Which shows that even when elements of weight matrix M are of order higher than 1/yN, 
density function of J A + M still converges to that of J A as iV — > oo. This is due to the column 
structure of M. In comparison, we show in panel 4 the density function of W = J A + M, where 
M is a constant matrix whose elements are randomly chosen from uniform distribution on [0,1]. 



In FiglJl we compare the eigenvalue spectra of W — JA + M and W = J A. In both 
cases J is drawn from a non-Gaussian ensemble. All spectra are generated by Monte-Carlo 
simulations. Since p(x,y) — p(\z\), where z — x + iy, it is sufficient to show the dependence 
of eigenvalue density function on radius \z\. We find these functions match quite well. In 
comparison, we replace M with a constant matrix which does not have the column structure, 
and find the density function is rather different. 

Next, we test our result on the two-component Gaussian network of pQ. Denote the weight 
matrix by W — J A + M and follow the brief discussion in the introduction. We fix / = 0.2 
and choose the variance of the Gaussian distributions for excitatory and inhibitory synaptic 
strengths to be and respectively. The only constraint on the mean strength matrix 
M is the 'balance' condition. As in the previous example and already pointed out in pQ, the 
bulk of the spectra of W with different M matrices are almost identical. From numerical 
simulations, it appears that when M has larger elements, there are more eigenvalues outside 
the support of the spectrum, see Eq. (|3.22p for the definition. To show this is a finite-size 
effect, we calculate the percentage of 'outliers', for different matrix size N. It shows in Figj3] 
that, for all values of M, as N — > oo, percentages of 'outliers' for W = JA + M approach to 
that of W = J A. 
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Figure 3: Percentage of eigenvalue outliers of two-component Gaussian network W decays as N 
increases. For / = 0.2, elements in fN columns of W have Gaussian distribution with variance 
Ty. and mean \x e . Other elements have Gaussian distribution with variance jr and mean in. 
From top to bottom, (ji e ,Hi) are (Ij-j)) (jf^~w) amd (°'°)' respectively. 

3 Synaptic strength of Gaussian distribution 

In this section, we calculate the eigenvalue density of multi-component Gaussian network. 
Assume there are m types of neurons in the network and synaptic strengths have different 
Gaussian distributions. From the previous section we know the density functions of J A + M 
and J A are identical when M has the column structure and satisfies the 'balance' condition, 
even without the additional constraint of pQ. Therefore, spectrum of this network is the 
density function p(x, y) of the following random matrix 

W = JA, (3.17) 

where J is drawn from the real Ginibre ensemble, i.e. V(x) = \x in Eq. (|2.4[) . and A is 
defined in Eq. (|2.1l) . The case to = 2 is solved in [1] with the method in [6]. For to > 2, we 
find the technique developed in [22, 23 is more convenient. Define an operator tr^r which, 
when acts on an N x N matrix A, gives trjv^4 = tr^^lA 2 . By Eq. (|2.15j> . the equation for the 
one particle irreducible (1PI) self-energy T,w is 

y _ ( Si s 2 \ _ i / o tv N g 2 -i N \ . , 

^ W ~{^ S 4 J-^ v tr^ 3 -A 2 )■ ^ 

Note that the above matrix has the structure outlined in Eq. (|2.16j) . The generalized Green's 
function Qw (defined in Eq. (|2.9l) ') is related to T,\y by the Schwinger-Dyson equation Eq. (|2.10[) . 

<*>-{% ZHi-Z i:g)l„- 



From Eqs. (|3.18p and (13.191) . we get the following equation for an unknown variable p 



E 



z\ 2 -paf 



1. (3.20) 
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Figure 4: Self-energy H\y is related to resolvent Qw by equation Eq. (|3.18p . 



Eq. (|3.20p has multi- number of solutions. The correct one for our problem is the one satisfying 
the boundary condition 



-1. 



(3.21) 



The boundary of spectrum is determined by the transition point p — |22j , which corresponds 
to the circle with radius \z\b, such that 



(3.22) 



The disk region defined by \z\ < \z\b is the support of the spectrum. Off the support, we 
always have p = 0. In the case m = 2, the formula in Eq. (|3.22| > gives the same result for 
spectrum boundary obtained in [1] by solving a saddle-point equation. From (I2.11j) . the 
Greens's function Gw{ z , z) is given by the following formula 



G w (z, z) 



Eva 
i=l 



fi 

\z\ 2 —pcr'f '■ 



z z < z 



\z\ A > z 



From Eg. (13.201) . when on the support of spectrum, we get 



_ i=l (\z\ 2 -p<r?) 2 



d\z\ 2 1 



Finally, by Eq. (|2.7p . we get the eigenvalue density of W 



p(x,y) = 



1 (\ z \2_9p_ _ \ V^m ha_ 

7T \ \^\ d\z\ 2 l^ii=\ (| z p-p 



0. 



z?< z 



kl 2 > 1*1! 



(3.23) 



(3.24) 



(3.25) 



Introduce the notation (a a ) = h a ii where a is a constant. From Eq. (|3.25p . we find the 
eigenvalue density at the centre and boundary of the spectrum 



P(0) 



(a 2 ), and p{\z\ B ) = 



1 (a 2 



(3.26) 



When m = 1, from Eqs. (13.201) and (|3 . 25[) we easily recover the well known result for 
Ginibre ensemble [12]. For m = 2, choosing the solution for quadratic equation Eq. (|3.20[) 
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Figure 5: Density p of eigenvalues as a function of radius in the complex plane |z|, for 
N = 400. The solid lines are the analytic results by Eq. (|3.25p and symbols are numerical 
simulations. The figure shows results for different sets of variances a 2 /N with fixed population 
/ = (0.1,0.2,0.3,0.4). 
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Figure 6: Density p of eigenvalues as a function of radius in the complex plane \z\, for N = 400. 
The solid lines are the analytic results by Eq. (|3.25p and symbols are numerical simulations. The 
figure shows results for different sets of population with fixed variances a 2 = (0.1, 0.2, 0.3, 0.4) /N. 
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satisfying condition Eq. (|3.21[) . then from Eq. (|3.25j) . we successfully recover the results ob- 
tained in pQ. For large m, it is hardly possible to have an analytic solution for Eq. ([3.20(1 . 
But it is very simple to find the numerical solution for this algebraic equation. It turns out 
that there is always only one solution on [—1,0], which is just what we need according to 
the boundary conditions. 

In FigEl we compare the density function in Eq. (|3.25j) with numerical simulations for 
synaptic strengths with different variances er's but the same /'s. In Fig|6j we let fs change 
but keep a's fixed. In both cases, we observe very good match between numeric data and 
analytical results. The only significant deviation happens near \z\ — 0. In fact, this deviation 
already appears when A = In and is shown due to finite-size effect [6]. 

4 Discussion 

In the first part of this paper, we show that modifying the mean strengths of synapses of 
a neural network does not change the density function of synaptic matrix even when there 
are several types of neurons and the strengths of their synaptic connections have correlated 
distributions. 

In Eq. (|2.4p . the ensemble of random matrix J is chosen to be O(N) invariant so that 
all elements of random matrix J have the same distribution. Differences between different 
types of neurons are introduced only by A and M. In fact, we can draw the synaptic matrix 
W from more general ensembles. As long as Eq. (|2.15p holds and M has the column block 
structure, eigenvalue spectra of W will not depend on M. 

We therefore prove that the density functions of large random matrices described by 
Eq. ([2.3|) - ([2.4|) are not changed by perturbations of the rank-1 matrix M. This type of 
random matrices are not of Wigner type or have iid elements as in jTTJ [TB] . 

It is its structure that makes M irrelevant to the eigenvalue density function. In reality, 
we may need to choose the mean value of synaptic connections to be of the same order 
of their fluctuations, i.e. (My) cx N^ 1 / 2 . But this is not necessary in our proof. If the 
'balance' condition is not imposed, the eigenvalue spectra will be identical to the 'balanced' 
case except the eigenvalues at zero will be shifted [TJ [T7] . 

In the second part of this paper we calculate the density function of random matrices 
of the form J A, where J belongs to Ginibre ensemble. These matrices describe random 
networks with multiple independent components. We find closed formulas for the eigenvalue 
density at both the centre and the boundary of the spectrum in terms of variances of synaptic 
strengths. 

When J is drawn from the ensemble in Eq. (|2.4p and A = In, we know by the Single- 
Ring Theorem [24l [25l [26] that the support of the eigenvalue spectrum is either a disk or 
an annulus. It will be interesting to find out whether or not the Single-Ring Theorem still 
holds when A is diagonal but not proportional to In ■ Eq. (|3.22[) shows when J has Gaussian 
distribution the support of the spectrum is always a disk of radius \z\b, but never an annulus. 
This indeed agrees with the Single- Ring Theorem. Clearly, to prove the Single- Ring Theorem 
for J A, where J belongs to the general ensemble defined in Eq. (|2.4l) . we need to take different 
approaches. Work on this topic is currently in process. 
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Appendix 



In this section we show that in the large N limit, due to the structure of A4, the Green's 
function Gw defined in Eg. (12.111) equals to Gj\ defined similarly for J A, i.e. 



G w =G./a + 0(7V-1). 



(A.27) 



Step 1. Define D = (z — M T )(z — M) — abA 2 , where a and b are scalars shown in Eg. (12. 161) . 
By Eg. (12.131) . Gw can be written as 



* (z-M)D- 1 



(A.28) 



Here we used the fact that the (12)-element of £ja is proportional to identity matrix and 
the following formula from linear algebra 



E F 
G H 



E~ 



E^FX^GE- 1 
-X^GE- 1 



-E^FX- 1 

x- 1 



(A.29) 



where X = H - GE^F. 

Step 2. Let I pq be a p x g-dimensional matrix with all elements equal to 1 and let 
Mi = fiN for i = 1, . . . , m. Then D^ 1 has the following m x m-block structure 



D- 1 = 



( ax + bulM-iMi 

&2iIm 2 Mi 0,2 + &22IW2M2 



b2m^M 2 M„ 



\ 



(A.30) 



+ b mm lM m M m ) 



where a, = l/(\z\ 2 - ab\ 2 ), b u - O^ 1 ), = 0(N- 2 ), and a, + b lt M t = OiN- 1 ), 
for i,j — 1, . . . ,m. This claim is proved by induction. First, notice that D has the same 
m x m-block structure as in (|A.30[) except all its parameters are of order 1. Using the fact 



(a + bl M Ai) 



1 



a a{bM + a) 



MM, 



(A.31) 



and (|A.29I) . we find that for m = 2, D^ 1 indeed has the properties described by (|A.30I) . 
Then assume the claim is true for m = n — 1. By straightforward calculation using (|A.29[) 
and (| A.31)) . we find the claim is also true for m = n. 

Step 3. Substituting (|A.30|) to (|A.28j) and using Eq. ([2~TT]) . we get (|A.27j) . This completes 
the proof. 
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